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We present a new code for solving the molecular and atomic excitation and radiation transfer problem in a molecular gas and predicting 
emergent spectra. This code works in arbitrary three dimensional geometry using unstructured Delaunay latices for the transport of 
photons. Various physical models can be used as input, ranging from analytical descriptions over tabulated models to SPH simulations. 
To generate the Delaunay grid we sample the input model randomly, but weigh the sample probability with the molecular density and 
other parameters, and thereby we obtain an average grid point separation that scales with the local opacity. Our code does photon very 
efficiently so that the slow convergence of opaque models becomes traceable. When convergence between the level populations, the 
radiation field, and the point separation has been obtained, the grid is ray-traced to produced images that can readily be compared to 
observations. Because of the high dynamic range in scales that can be resolved using this type of grid, our code is particularly well 
suited for modeling of ALMA data. Our code can furthermore deal with overlapping lines of multiple molecular and atomic species. 
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1. Introduction 

The ability to predict the radiation signature associated with a 
physical model is crucial when interpreting astronomical ob- 
servations. In the general case, this can be achieved by solv- 
ing the equation of radiation transport and, in the special case 
of non-local thermodynamical equilibrium (non-LTE) radiation, 
the equilibrium quantum state distribution. This is a strongly 
coupled problem that can only be solved by numerical analy- 
sis except for the limiting cases where either the material is in 
LTE or in the optically thin regime where there is little interac- 
tion between matter and radiation. The computational demand 
increases rapidly with complexity of the model if such approxi- 
mations cannot be applied. Due to the low density conditions of- 
ten found in inter- and circumstellar material, LTE is oftentimes 
not a valid approximation and therefore, if the medium turns op- 
tically thick, numerical methods must be used to calculate the 
molecular excitation in such environments. 

Ever since the first molecules where discovered in 



space (CH + and OH: |Douglas & Herzberg|194T||Weinreb|1963| ), 
line emission has been an important astrophysical tracer of the 
physics and chemistry. The chemistry can be traced from the 
abundance of various molecular species and the physical envi- 
ronment, such as temperature and density, from the excitation 
of the lines and kinematical information is embedded in the 
line profiles. In star forming regions, molecules are also trac- 
ers of shocks and can be used to investigate the external radia- 
tion environments ( [van Dishoeck et aL| |2009). As an example, 
the inside-out collapse model of low-mass star formation ( |Shu| 
|1977| ), which is often used by many authors, has been confirmed 
for a number of objects by measurements of vari ous molecular 
lines (e.g., |Walker et al.||1986| |Zhou et al.||1993| ). However, in 
order to derive reliable molecular abundances, it is important to 
understand the spatial distribution of the gas, and while the gas 
distribution in the plane of the sky is directly measurable, given 



sufficient instrument resolution, the distribution along the line of 
sight is not easily measured. In fact, deriving the three dimen- 
sional distribution based on a two dimensional image forms an 
inverse problem similar to the problems encountered in medical 
tomography and seismology. In the astrophysical case, we deal 
with such problems by constructing a source model based on our 
theoretical understanding of the object in question, and then cal- 
culate the line signature of this model which is then compared to 
the observations and, based on the result, the model is improved 
to give a better fit. 

A number of numerical codes and algorithms that focus on 
predicting molecular emission in the far-infrared and millimeter 
wavelength regimes hav e been made available to the community 
in the last decades (e. g. , |B ernes 1 1 97 9 [ | Juvelaj 1 99l\ [Ho gerheij de 
| & van der T^|20U0l |Schoier||2000[ |Pavlyuchenkov & Shustov 
2004J). As computers have become faster, codes have evolved 
from working predominantly in a spherical symmetry (ID) to 
being able to solve problems in a cylindrical symmetry (2+ ID). 
This progress has mainly been driven by the development in tele- 
scope facilities, providing observers with ever increasing spatial 
resolution. This paper describes a new non-LTE radiation trans- 
fer code for (sub-)millimeter and far-infrared continuum and 
spectral line radiation that works in arbitrary geometries (3D). 
The code has been designed to be fast, reliable, and easy to 
use, with particular emphasis on being able to solve models in 
very high resolution for the Atacama Large Millimeter Array 
(ALMA) and models of molecules with a complex level config- 
uration (e.g., H2O and CH3OH), relevant for observations with 
the Herschel Space Observatory. The code is called LIME (Line 
Modeling Engine). 

LIME derive s from the RATRAN code (Hogerheij de & 
van der Tak 2000) and, although rewritten from scratch, it shares 
much of the code base and the solution method. RATRAN is 
a Cartesian grid based Monte Carlo code which includes el- 
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ements of accelerated Lambda iteration (Rybicki & Hummer 
|1991| ) and it exists in both a one and two dimensional version. 
The main difference between the two codes is the photon prop- 
agation method (see Sect. |3.1| ) which dramatically improves the 
calculation time, making 3D models feasible to solve, and al- 
lows for a very flexible model input method. Our photon trans- 
port method is a line extension of the approach of |Ritzerveld 
|& Icke| ([2006 ) for their continuum radiation transfer algorithm 
simplex. This algorithm was originally developed in to provide 
a fast solution to the radiation transfer problem in three dimen- 
sions for all opacity regimes, in particular for hydrodynamical 
cosmology simulations, where approximate radiation transfer is 
often used to speed up temperature calculations. The simplex al- 
gorithm has been shown to work well for clumpy and inhomoge- 
nous media with a large density contrast. LIME has been written 
mainly for the purpose of solving models of disks and envelopes 
around young stellar objects (YSOs), but due to the 3D nature 
of the code, it can, however, also solve models of more complex 
structures such as outflows, molecular clouds and even clusters 
of YSOs and their environments. The main limitation to this is 
the availability of input models. It is to be expected that mod- 
els provided by (magneto-) hydrodynamical simulations will be- 
come more important in the near-future, and it has therefore been 
important for us to provide an easy way of interfacing the out- 
put of a dynamical simulation with this code. SPH simulations 
are in particular well suited as input for LIME, since the physi- 
cal quantities is already described by particles in very much the 
same manner as LIME uses it. 

A number of new features have been added to LIME. As it 
is expected that ALMA will detect potentially hundreds of lines 
per setting and many of these lines may be blended, we have 
made it possible to solve the radiation transfer for an unlim- 
ited set of molecular species simultaneously and dealing with 
the cross-excitation correctly. This feature obviously comes at 
a cost of increased calculation time, of which the severeness de- 
pends on the number of coupled lines. Another new feature is the 
extended output capabilities, where the user can ask for spatial 
intensity and opacity distributions as well as the grid information 
be written out for visualization purposes. 

The outline of this paper is as follows: Section [2] describes 
the physical problem of radiation transport in a molecular 
medium, Sect. [3] provides an extensive description of the LIME 
code, Sect. [4] brings a few examples and validation of the LIME 
results, and finally, conclusions are given in Sect. [5] 



2. The physical problem 

The spectral intensity I v of radiation propagating through a 
medium with emission and absorption coefficients j v and a v is 
described by the equation 



dly _ . 

ds ^ v 



(1) 



where we neglect the scattering term which is of little relevance 
in the sub-millimeter regime. The coefficients j v and a v are re- 
lated to the Einstein coefficients A\ 2 , B\ 2 , and B 2 \, for transitions 
between any two adjacent levels with an energy of hv = E\- E 2 , 
by the equations 



hv 

Jv;gas = ^n 2 Ai2(f>(v) 
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Vv, g as = —(niB u -n 2 B 2 i)(p(y), 



(2) 
(3) 



for the line radiation and 

jv;dust ~ ~®v;dustBv(T dust) 
&v;dust = ^vPdust-> 



(4) 
(5) 



for the thermal dust emission, (p is a function of frequency which 
contain the contribution of several line broadening mechanisms 
and h is Plancks constant. The dominating mechanism here is 
the doppler broadening due to local turbulence. B v is the Planck 
function for a given dust temperature and k and p are the dust 
opacity and dust mass density respectively. The dust opacities 
depend on the type of dust and various tab ulated descriptions can 
be found throughout the literature (e.g., Ossenkopf & Henning 
1994]) The contribution from the gas and the dust components 
add up to form the total emission and absorption coefficients. 

At any given position, the local mean radiation field can be 
obtained by solving Eq. [T] and integrating the intensity over all 
solid angles 
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The gas can be excited either through absorption of a photon 
or through collisions with other gas molecules. Collision rates 
between two energy levels i and / depend on the local kinetic gas 
temperature 



Cu = ^c a exp(- — ^ (Ei - Ed |, 

gas 



(7) 



where g t is the statistical weight of the /'th level and cu is 
the molecule dependent rate coefficients, which, in general, are 
functions of temperature as well, ks is Boltzmanns constant. The 
molecular excitation thus depends on Ju (where il replaces v to 
denote the frequency associated with transition between level i 
and /) and Cu and the fractional population of the /'th level ni is 
given by 



Hi 



^i>in k Au + Y. m ni(BuJu + C u ) 
^i<Au + ^m(BuJn + Cu) 



(8) 



assuming statistical equilibrium. Because Ju implicitly depend 
on Hi through j v and a v , Eq. [6] and [5] forms a recursive problem 
and must be solved iteratively. The collisional rates cu in Eq. [7] 
are not provided as part of the LIME package and must be ob- 
tained separately. LIME, however, can read the collisional data 
files provided by the LAMDA databas^] 

The whole problem is further complicated if a systematic 
velocity field is present. In this case photons will be Doppler 
shifted and thereby contribute to different transitions if the tran- 
sition energies are spaced closely enough. While this is not a 
problem for molecules with sufficiently widely spaced transi- 
tions, such as CO, more complex molecules such as CH 3 OH 
have a much richer level sub-structure where photons can eas- 
ily contribute to several transitions given a systematic velocity 
field or even by random velocities in the gas. Indeed, a photon 
originating from a transition in one molecule may get Doppler 
shifted and excite a transition in another molecule, and there- 
fore we track photons by their frequency, not relative velocity, 
for each transition in question, summing up the contribution to 
J v from all other sufficiently close transition. 
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3. Solution method 

3.1. Computational grid 

In radiation transfer codes, the source model is typically laid out 
onto either a grid of rectangular cells or an (r, #)-grid. Model 
properties such as density, temperature, molecular abundance, 
as well as the radiation field J v and the level populations n t are 
taken to be constant over each cell. Photon packages are traced 
in random directions from random points of absorption within a 
cell and Eq. [T]is integrated along the paths. Summing over all 
photons, the mean radiation field J v is obtained for the cell. 

In our implementation, the source model is not mapped 
onto a regular grid of cells. Instead we use a random set of 
points which represents the local environment (density, tem- 
perature, populations, etc.). The points are distributed in three- 
dimensional space, and thus we are able to use three-dimensional 
source models whereas in, for instance RATRAN, even though 
photons propagate in three-dimensional space, the source model 
needs to be rotational symmetric around the second axis and 
mirror symmetric around the first axis. Our points are placed 
randomly throughout the entire computational domain, however 
with a probability that is weighted by a source model depen- 
dent function. This approach is similar to the one described 
by |Ritzerveld & Icke| ( [20 06). We choose to use the molecular 
density profile of the source model as a probability distribution 
for the grid points, so that we end up with a point distribution 
that has particularly interesting properties from a radiation trans- 
fer point of view, namely, that the average distance from a point 
to its neighbors becomes inversely proportional to the density 
and proportional to the local mean free path, since the mean free 
path / = (a v p)~ l . 

If we at first consider continuum radiation transfer only, 
where the absorption coefficient is independent of the radiation 
field, a v is constant and the mean free path depends on density 
alone. We can thus obtain a grid point distribution where the ex- 
pectation value of the neighbor point separation equals the local 
mean free path by adjusting the number of grid points in accor- 
dance to the dust opacity k v . 

Figure [T] shows that the average neighbor distance is propor- 
tional to the density over two orders of magnitude. The arbitrary 
offset between the mean free path graph and the point separation 
graph scales with the total number of grid points. The point sep- 
arations are seen to deviate from proportionality at radii smaller 
than 10 12 m (« 7 AU). This is because we impose the constraint 
that a certain number of points should be present a the largest 
scales and thus it would take an unfeasible number of points 
to follow proportionality with the mean free path down to 0.1 
AU. The point distribution in Fig. [T] consists of 40000 points. 
LIME offers several different options for sampling the density, 
depending on the geometry of the input model, such as uniform 
sampling in a rectangular box or logarithmic radial sampling in 
a sphere. It is also possible to use a fixed set of points, for in- 
stance the cell centers of a tabulated input model or the particles 
of an SPH simulation. In this case, the resulting point distribu- 
tion may not, of course, scale with the mean free path. In other 
cases it may not be desirable to grid according to the molecular 
density. If, for instance, the user is interested only in the very 
high /-lines of a molecule, it can be useful to grid according to 
the temperature distribution, or if the model involves shocks or 
outflows, the velocity gradients. In that case however, the grid 
does not describe the opacity well, but rather ensures that cer- 
tain spatial region are well sampled. It is possible to logically 
OR point distributions, so that a grid may be based on both the 
temperature and the density, in which case the opacity is also 




Fig. 1. The dashed line shows the local mean free path in the 
continuum for a model with a density profile p oc r" L5 and a con- 
stant absorption coefficient. The full line is the average neighbor 
distance for a random density weighted point distribution. 



sampled. In fact, point distributions for each individual opacity 
source (e.g., gas and dust) should always be OR'ed to form a 
single grid. 

If we now consider line radiation transfer, the local ab- 
sorption coefficient a v depends on the current level population 
through Eq. [3] The local mean free path thus changes with the 
radiation field and no longer scales simply with the density, but 
rather with the population dependent opacity. Furthermore, the 
opacity varies across the line, so photons at different frequen- 
cies will have a different mean free path, and therefore the local 
photon mean free path is not a well defined property. Despite of 
this, gridding according to density turns out to be the best way 
to form the grid for line radiation transfer too, simply because 
it results in a grid that describes the spatial distribution of the 
molecules very well. 

When the point distribution has been obtained, the points are 
connected by Delau nay triangulat ion, using the public available 
QHull library (Barb er et al.|1996| ). In 2D, the Delaunay triangu- 
lation is constructed by connecting any three points that defines 
an empty circumcircle, meaning that no other points can lie in- 
side the circle defined by a Delaunay triangle. The definition is 
valid for all higher dimensions and in 3D the Delaunay trian- 
gulation forms tetrahedra out of four vertices. Figure [2] shows 
a random point distribution (in 2D) in the left panel with its 
Delaunay triangulation shown in the center panel. Also shown, 
to the right, in Fig.|2]is the corresponding Voronoi diagram. The 
Voronoi diagram is the topological dual of the Delaunay triangu- 
lation and one can be constructed from the other. For a discrete 
set of generating points P, a Voronoi cell is defined as the region 
surrounding the site s e P which contain points that lie closer 
to s than to any other generating sites in P. The physical proper- 
ties of the grid points (density, temperature, excitation, etc.) thus 
represent the entire Voronoi cell associated with that point and 
therefore the Voronoi cell can be considered similar to the cells 
in traditional codes. As a matter of fact, if we chose a regularly 
sampled and unweighted grid point distribution, the Voronoi di- 
agram will describe an ordinary Cartesian mesh. One beneficial 
property of the Voronoi grid is that it automatically conserves 
the mass when mapping a density function. This is difficult to 
do on a Cartesian grid because the mass centroids of the cells 
are dependent on the cell orientation with respect to the underly- 



4 



C. Brinch and M. R. Hogerheijde: LIME - a flexible, non-LTE line excitation and radiation transfer method 




Fig. 2. The leftmost panel shows a random point distribution. In the middle panel, the points have been Delaunay triangulated. The 
rightmost panel shows the corresponding Voronoi tessellation. 



ing model. Furthermore, because of the random orientation and 
shape of its cells, Voronoi grids do not suffer from the aliasing 
effects which are inherent to regular Cartesian grids. 

When the grid points have been distributed throughout the 
model domain, it is inevitable, due to the stochastic nature of the 
sampling method, that some points end up much closer than the 
local separation expectation value and some will be much fur- 
ther apart. This results in a Delaunay triangulation that is very 
irregular with some triangles being very long and narrow. This 
irregularity can be remedied by applied what is known as Lloyds 
algorithm ( |Lloyd|19 82; Springel 2010 ), which iteratively moves 
a grid point slightly toward its Voronoi cells center of mass. In 
our implementation, each grid point is moved slightly away from 
its nearest neighbor for a preset number of iterations. The effect 
is illustrated in Fig. [3] for a random 2D sampling of a Gaussian 
density profile. The top panels shows the initial unsmoothed dis- 
tribution while the lower panels showed the triangulated point 
distribution after the smoothing algorithm has been applied. The 
plots in the right column show the neighbor distances as a func- 
tion of radius. In the smoothed grid, the distances are much less 
scattered and follows the Gaussian profile (shown as the light 
colored full curve) much more accurately than in the top panel. 
The smoothing strategy should not be exaggerated, i.e., moving 
the points too slowly and iterating trough too many stpdf, be- 
cause the algorithm will then act as an annealing process and it 
will result in a perfectly regular grid where all variations in the 
point distribution due to the underlying density field is smeared 
out. By doing it right, however, a smooth grid can be obtained 
while the underlying density structure is still preserved in the 
grid. We have found empirically that by using 25 iterations and 
moving the closest neighbors about 10% of the distance away 
from each other results in a sufficiently smooth grid that pre- 
serves the underlying physical structure well. 

During gridding of our source model, we also distribute 
a number of points randomly on the surface of a sphere sur- 
rounding our model. These points are also Delaunay triangu- 
lated and connected to the model grid points, but they do not 
represent anything except the surface of our computational do- 
main. Whenever a photon reaches one of these sink points, it is 
considered to have escaped the model. 
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Fig. 3. An unsmoothed (top) and smoothed (bottom) Delaunay 
grid based on a Gaussian density distribution. The right hand 
column shows the neighbor point separations as function of ra- 
dius. The yellow lines are aids to the eye. They show a Gaussian 
distribution that describes the density distribution. 



3.2. Photon propagation 

The photon transport itself goes along Delaunay lines only, from 
one point to another, which makes integration of Eq. [T] partic- 
ularly simple and very fast. In the three-dimensional Delaunay 
triangulation, the expectation value for the numb er of lines at- 
tached to a grid point is approximately 16 (Ritzerveld & Icke 
2006 ) and the spatial sampling of J v is thus limited to this num- 
ber of directions. However, we still need to trace a number of 
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photons along each Delaunay line, not only in order to sample 
the frequency band properly, but also because we cannot con- 
serve momentum stringently with a single photon on this grid. 
In principle, a photon passing a grid point from a certain direc- 
tion should continue to travel in the exact same direction. This 
is in general not possible due to the random orientation of the 
Delaunay lines, so instead we choose one of the two outgoing 
Delaunay lines (H and £2) that make the smallest angle with the 
original direction of the photon. The outgoing line is picked at 
random, but weighted by the ratio of the two angles, 



1.0 



p{ll) = (Z1/Z2WD, 



(9) 



where L\ < Z 2 . The same procedure is used at all subsequent 
grid points (using the original momentum vector to determine 
the outgoing direction) until the photons escape the model. By 
sending a number of photons along each initial Delaunay line, 
we thus probe, not a single line of sight, but rather a cone, while 
still conserving momentum on average. An example of the pho- 
ton propagation is shown in Fig. [4] for a single point and a sin- 
gle direction. Because of the relative low number of photons 
needed to probe the spatial directions, we can allow ourself to 
increase the number of photons used to sample different frequen- 
cies, while we still maintain a low (initial) number of photons per 
grid point. The inset in Fig. [4] shows the distribution of the lo- 
cation where the photons reach the surface of the grid. This dis- 
tribution is reasonably well described by a Gaussian distribution 
around the intersection of the original momentum vector and the 
surface. The number of initial photons is a user-defined setting, 
but as a default value, we use five times the number of neighbor 
points, so that each neighbor is initially probed by five photons. 
These photons are distributed evenly across a frequency range 
of +3<x with respect to the line center so that the median pho- 
ton coincides with the local rest frequency, cr is determined by 
the local turbulent velocity dispersion through the user-defined 
Doppler b-parameter. 

Any given grid point will see more Delaunay connections 
coming from high density regions than from low density regions, 
simply because the grid point density is higher in high density re- 
gions. Because of this inhomogeneity in the angular distribution 
of Delaunay connection, care must be taken when averaging the 
radiation field using Eq. [6j In our implementation, this equation 
reduces to a discrete sum 



1 ^£ 



(10) 



i=0 



where N is the number of Delaunay neighbors. o)i is a weight 
that is proportional to the solid angle represented by the fth 
Delaunay line. This angle corresponds strictly to a surface area 
on a unit sphere, but we use the area of the Voronoi facet that cor- 
responds to the Delaunay line as a good approximation (within 
10%). 

Figure [5] shows a comparison of J v between LIME and 
RATRAN. The input model is a thin flat disk with a density pro- 
file oc r~ l . The radius is 500 AU and the height is 50 AU. The 
disk is placed in an ambient low density field, n=10 4 cm" 3 . For 
LIME the model is sampled by 8000 points, whereas RATRAN 
uses 400 grid cells. The radiation field of first three levels are 
shown in panel a)-c) in Fig. [5] The LIME points are shown as 
black points and the RATRAN points in yellow. The points for 
each transition makes up two distinct populations, an almost 
horizontal branch and a scattered population below. The tight 
horizontal distribution of points are the ones that lie inside the 
disk. These points are extremely well matched between the two 
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Fig. 4. The propagation of photons from one point in one partic- 
ular initial direction. As the photons step through the grid they 
choose the following step by weighing the probability with the 
inverse angle the direction makes with the initial direction, i.e., 
a large angle produces a small probability of going in that direc- 
tion. The inset shows the distribution of photon arriving at the 
surface. A Gaussian distribution is overplotted for comparison. 



codes. The scattered population of points are the ones that fall 
outside of the disk radius and these are also well matched. The 
LIME points in the ambient low density region scatters a bit more 
than the corresponding RATRAN points do. This is not a dilution 
of the radiation field due to poor spatial sampling or erroneous 
photon propagation on the Delaunay grid, but simply because 
the LIME grid has a much higher resolution than the RATRAN 
grid. The proof of this can be seen in panel d) in Fig. [5] where 
a similar comparison of the / = 3-2 transition between LIME 
and RATRAN has been made, but with a LIME grid which is not 
density weighted at all, which means that all spatial regions are 
equally well sampled. This distribution is indistinguishable from 
the one in panel c). 

In LIME we integrate Eq. [T] all the way to the edge of the 
model at which point we add the external contribution, in most 
cases the cosmic microwave background. This is know as the 
method of long characteristics. Another approach exists, known 
as the method of shor t characteristics (Kun asz & Auer| [l988; 
|Auer & Paletou|p"994| ), where photons are only traced to the 
neighbor grid cell at which point an interpolated radiation field is 
added. Short characteristics could be implemented in LIME with 
a possible gain in calculation time. In fact, when we average the 
radiation field by splitting the photon packages along two outgo- 
ing Delaunay lines, it is reminiscent of the short characteristics 
method. 

The level populations n\ in a grid point, as defined by Eq. [8] 
represent the molecular excitation in the entire Voronoi cell as- 
sociated with the point. J v is then the radiation field seen by the 
grid point. Following the method of Rybicki & Hummer (1991), 
when working with Cartesian grids, the radiation transfer can be 
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Fig. 5. A comparison of J v between LIME and RATRAN. The 
small black dots are the values from the LIME code and the yel- 
low dots are the cell values from RATRAN. The blue dots in panel 
d) are also from a LIME model, but where all the grid points have 
been distributed randomly over the model domain (with no den- 
sity weighting). 



formulated iterative as 

4 +1 =A[Sy(4)l (11) 

where S v = j v /a v is called the source function and k refers to 
current iteration. A is a matrix in which each element represents 
the radiation coupling between each pair of grid points. Thus the 
diagonal represents the radiative contribution from the cell itself 
while the off-diagonals are the coupling between the cell and 
its neighbors. If the entries in the matrix are organized so that 
cells that are immediate neighbors lie close to the diagonal, it is 
possible to split A in the following way 



K xt [S v (J k v )]+A local [S v (J k v +1 )l 



(12) 



where Ai oca i is a di- or triagonal matrix (depending on the imple- 
mentation) representing the radiative interaction from the local 
neighborhood. This will entirely dominate over the external con- 
tribution in high opacity problems . Because the A\ oca i matrix is 
much easier to invert, this approach is much faster than solving 
Eq.fTTI directly and this is in essence what is kn own as acceler- 
ated A-iteration (ALI) ( [Rybicki & Hummer|199T1 ). 

Because we cannot order the matrix elements of A in a mean- 
ingful way, we do not actually construct the A operator, but 
rather split J v j tai = J v, local + Jv,ext> separating the radiation field 
in the parts that comes from the inside and the outside of the 
Voronoi cell for which we are solving the populations (Fig. [6]). 
This split results in a solution to Eq.[T]that looks like the follow- 
ing, 



J v = I v e- J + 5 V (1 - O. 



(13) 



Because S v is a function of j v and a v , which are evaluated using 
the populations in the previous iteration, this method is similar 
to ALI. r is the opacity of the cell itself, so that if ds denotes the 
distance from the grid point to the Voronoi face r y = a v ds. In 
the case where the face of the Voronoi cell is located exactly one 
mean free path away from the grid point, 
-l 




Fig. 6. The radiation field seen by a grid point i is split into a lo- 
cal part J i oca i that originates from within the /'th Voronoi cell and 
an external part J ext that originates from outside the /'th Voronoi 
cell. The size of the Voronoi cells, i.e., the region in which the ra- 
diation is considered "local", scales with the inverse of the den- 
sity. 

This is exactly the condition we are striving to obtain by gridding 
according to the density as described in Sect. |3.1| Of course, the 
mean free path is frequency dependent, as discussed above, and 
therefore we can only achieve a condition where r scales with 
inverse of the density. The populations ni can now be obtained 
from the current J v by solving the linear set of equations given 
byEq.[8| 

3.3. Convergence 

One of the key problems when obtaining an equilibrium solu- 
tion iteratively is to decide when the solution has converged. 
There are many ways in which we can check for convergence, 
but it is necessary to find a method that applies to all models 
and molecules. It is important that the code does not quit pre- 
maturely before the model has actually converged, but on the 
other hand, it should not continue indefinitely because of minor 
random fluctuations in a single grid point somewhere. Given a 
sufficiently large number of grid points, the random nature of 
the code implies that some points will always deviate from the 
equilibrium populations. In LIME we therefore use a statistical 
criterion, and let the user decide how many iterations he or she 
will let the code run. We have found empirically that 15-20 iter- 
ations a good value for most models and the default is thus set to 
20. This can however easily be adjusted by the user according to 
preference and the problem at hand. 

Figure [7] shows, as an example, the convergence history of 
a single grid point in a test run. Panel a) shows the signal-to- 
noise ratio of the current iteration, where the signal is the current 
population and the noise is the standard deviation of the previous 
five iterations (shown in panel b), 



S/N = 



(15) 



r y = a v (a v p) oc 1 1 p. 



(14) 



for each level. The signal-to-noise ratio is seen to increase and 
then level out after about 10 iterations. It does however level out 
at different values for the different levels which makes it diffi- 
cult to fix a certain value to reach. Also, the signal-to-noise ratio 
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Fig. 7. This figure shows the convergence history of a single grid point in a LIME run. Different colors refers to the quantum states 
0-7. Panel a) shows the signal-to-noise ratio of the current iteration. Panel b) shows the standard deviation, calculated from the 
previous five iterations. Panel c) shows the average population over the previous five iterations and panel d) shows the derivative of 
the populations. 



fluctuates a lot (notice the log scale) which again makes it hard 
to decide whether or not the solution is stable. In panel c) we 
can see the populations averaged over the previous 5 iterations, 
and from this plot it is quite obvious that a stable solution (for 
this grid point) has been reached after the 12'th iteration. The 
derivative of this curve is shown in panel d). 

Because the signal-to-noise fluctuates randomly from one 
iteration to the next and the range in signal-to-noise values is 
large, we have chosen to consider the distribution throughout 
the entire model. The median value of this distribution tells us 
about how well converged the model is in general. We also con- 
sider the minimum value for the signal-to-noise for all levels 
and grid points. Figure [8] shows the signal-to-noise distributions 
for the energy levels 0-6 individually and for all levels with a 
fractional population higher than 10" 12 . We use this cut-off, be- 
cause levels which are less populated only add unwanted noise. 
The differently colored histograms show the distributions with 
increasing iteration number. On top of the distributions, the cor- 
responding median values are marked. All medians are seen to 
increase with increasing iterations with the median of the least 
converged level 2 ending up at a signal-to-noise of 200 already 
at iteration 16. Still, the lowest signal-to-noise value of the entire 
model is as low as 20, which means that for that particular level 
we make an error of at most 5%. We find this acceptable and 



therefore stop the calculation at this point. However, for better 
confidence, more photons and iterations can be used at the cost 
of calculation time. 

Another test which the user can perform in order to evaluate 
the convergence, is to plot the spatial distribution of the 10% 
least converged points (or, for instance, all points with a signal- 
to-noise less than 100) and make a statistical comparison with 
the global point distribution to see if the least converged points 
in any way are associated with particular regions of the model or 
if the least converged points simply are a random subset to the 
entire grid. 

One particular situation that the user should be aware of 
is the very highly opaque regime. In this regime the radiation 
transfer problem becomes similar to the diffusion problem and 
this can cause a very slow drift of the populations toward the 
correct solution. This drift may be so slow that the populations 
seem converged and this problem is inherent to radiation transfer 
codes. 

3.4. Ray -tracing 

After convergence has been reached, the code ray-traces lines 
of sight through the model in order to obtain an image cube of 
the radiation that escapes from the surface. The user provides 
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Fig. 8. Distributions of signal-to-noise ratios in a test run for energy levels - 6 as well as for all significantly populated levels. The 
color coding refers to different iteration number. The thin vertical lines mark the median of the distributions. 
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the information on the source distance, source velocity, source 
orientation, and image resolution and units. The orientation pa- 
rameters are slightly more complicated than for 2D codes, where 
an inclination and position angle are enough to set the source ori- 
entation. In LIME we also have a source rotation which allows 
us to view a 3D model from any direction, using the matrix 



Red - 



cos((p) -sin((p) 

sin(6)sin((p) cos{6) sin(6)cos((j)) 
K cos(6)sin((p) -sin(0) cos(6)cos((p) 



(16) 



where 6 is the traditional inclination (0: face-on, n/2: edge-on) 
and (p is the azimuthal rotation. Rotation in the image plane is 
done afterwards, simply by rotating the image cube. 

For the ray tracing we let the photons move in straight lines, 
rather than jumping from grid point to grid point. In this part of 
the code we therefore do not make use of the Delaunay trian- 
gulation but rather the Voronoi diagram. The entire volume of 
the Voronoi cell is represented by the populations of the corre- 
sponding grid point and so integration of Eq.[T]becomes a matter 
of stepping through the source model and figuring out in which 
Voronoi cell the photon is. From an algorithmic point of view, 
this comes down to a simple sorting, not of thousands of cells, 
but only of, on average, the 16 neighboring cells. The step size 
is chosen as a fraction of the cell size in order to avoid acciden- 
tally missing a cell by stepping over it. Also, we need to sample 
variations in the velocity field across the cell in order to get a 
smooth spectrum. This is a very fast process compared to mov- 
ing through a regular grid, but not as fast as moving along the 
Delaunay lines, which is why this transport method is not used 
when determining the level populations. 

Apart from the images, which are output in standard FITS 
format, LIME can output a number of model diagnostics that are 
also obtained during the ray-tracing. This includes the opacity 
and intensity per Voronoi cell as seen along the line of sight. This 
information can be used to plot the r = 3 surface and identify 
the origin of the radiation for various transitions and tracers. It 
is also possible to dump the grid to a text file together with the 
populations and physical quantities so that the grid structure and 
excitation temperature as well as the input model itself can be 
explored visually. 

3.5. Benchmark: 1D collapse model 

A standard problem in radiation transfer benchmarking is the 
spherical Shu-collapse model ( Sh u|1977] ). This model describes 
a gaseous envelope that undergoes a gravitational inside-out col- 
lapse to form a protostar. In this example we use the particular 
fo rmulation of the problem fou nd in the code comparison project 
by |van Zadelhoff et aL| ( 2002 ) so that we may compare our re- 
sult to those of the codes that participated in that projecQ The 
molecule in consideration is HCO + . Eight different codes were 
tested against each other, both in terms of performance and their 
solutions to the problem. In the present test we will only com- 
pare the solutions because hardware and compilers have changed 
too much in the last ten years for a performance comparison to 
make sense. 

The model is already provided as a logarithmically spaced 
table describing 50 grid cells, shown in Fig. 9] Since LIME does 
not use regular grids, we interpolate this table and sample it ran- 
domly. For the comparison of solutions we have made a plot 
of the fractional level populations after convergence has been 
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Fig. 9. The model setup for the ID collapse benchmark test. 
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Fig. 10. Comparison of solutions to a ID collapse model be- 
tween different codes. The blue, red, and green curves are dif- 
ferent solutions from the RADTRANS code comparison project. 
The yellow curve is the RATRAN solution and the black line is 
the solution from LIME. 



reache d. There are two problems , denoted 2A and 2B in the pa- 
per by van Zadelhof f et al. ( 2002 ), corresponding to an optically 
thin and thick case. The fractional abundance of HCO + is 10" 9 
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and 10 respectively. We show the result of both tests in Fig. 
In this figure the fractional population of the ground state anc 
the first four excited levels are shown. All eight codes partici- 
pating in the RADTRANS comparison project agree very well. 
We have included only three of the eight solutions here (Doty, 
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Dullemond, and Juvela) as well as the solution of RATRAN in its 
present version. The black fluctuating line is the LIME solution 
and is seen to agree very well with the established codes. In our 
3D code, two points may be at the same distance from the center, 
but separated by a large angle and they may not be radiatively 
connected, especially not in high opacity models. These grid 
points may not see exactly the same radiation field because of 
the random nature of the grid and photon transport and therefore 
the level populations may not necessarily be exactly the same. 
There are no fluctuations in the results of the ID codes because 
only a single cell with a single solution exist at a given radius. 
If we take the RATRAN solution to be the expected solution we 
can calculate reduced^ 2 values, 



N 



1 y (n(r)uME ~ n(r) RATRAN ) 2 



(17) 



for each of the 6 levels of the LIME solution which gives us 
*io...5 ~ i 1 - 08 ' U1 ' L02 ' U7 ' L00 ' L0 °} for the optically thin 
case and x] =Q 5 ~ {1.06,3.12, 1.03,2.08, 1.28, 1.15} for the opti- 
cally thick case. The comparison is slightly worse for the opti- 
cally thick case, but here we also see greater variation between 
the established codes. 



4. Example 

The example we present is a typical 2D hydrostatic protoplane- 
tary disk model with a cold and dense mid-plane. Such models 
are numerous ly found i n the literature (e.g.,|Ch iang & Goldreich 
[T9971 [DuTIemond & Domimk|[2004l |Robitaille et al.|2006) but 
here we use a simple analytic toy model for illustrative purposes. 
The density structure is given by 



n(r,z) = n (r/r y L5 e- (z/h) \ 
where 



2Tk B r 3 



(18) 



(19) 



We consider HCO + , H2O, and CH3OH gas at a fractional abun- 
dance of 2xl0" 9 with respect to the H2 density. With these three 
species we illustrate the possibility of utilizing a large dynamic 
range in scales (HCO + ), very opaque models (H2O), and multi- 
ple overlapping lines (CH 3 OH). 

The temperature is given by a power-law 



T(r) = T (r/r ) 



-0.5 



(20) 



In a more realistic model the temperature would be calculated 
self-consistently based on the radiation properties of the central 
source and the temperature would drop toward the mid-plane be- 
cause this region would be shielded from stellar radiation by the 
upper layers of the disk. In our example we mimic this effect by 
lowering the temperature in a wedge shaped region around the 
mid-plane to 20 K. By letting water freeze out at temperatures 
below 90 K, we can simulate a com plex abundance structure of- 
ten used in protoplanetary disks (Jonkh eTd et al.|[2007 Woitke 
|et al.||2009| ). The disk extends to 500 AU and the values for n 
and T are 10 8 cm" 3 and 90 K at the radius of 100 AU. In ad- 
dition we have added a 2 AU wide gap around a radius of 5 
AU from the center. The disk is in Keplerian rotation. Figure 1 1 
shows the density, temperature, and H 2 density of our disk. Of 
other parameters that describes the disk are the turbulent velocity 



dispersion set to 150 ms" 1 , stellar mass of 1 M©, and a gas-to- 
dust ratio of 100. We use thin mantled grains with 10 7 years of 
coagulation and the resulting disk mass is 0.02 M©. 

To break the azimuthal symmetry and make the model fully 
3D, we have placed a protoplanetary condensation in the gap. 
The protoplane t has the same qualitative properties as the one 
described in Naray anan et al. ( 2006| ). The protoplanet has been 
modeled by placing a spherical distribution of grid points at the 
desired spot and giving the grid points an H2 number density of 
2 x 10 15 cm -3 , which, given a radius of 1000 Jupiter radii, results 
in a mass of about 1.4 Mj. The temperature of the condensation 
is kept at 150 K. 

Using this setup we have first made an edge-on view of 
the grid, which can be seen in panel d) of Fig. 11 The grid 



points (and their connections) are color coded according to den- 
sity, where blue is lower density and red is higher density. The 
number of grid points in this simulation is 10 4 for the disk and 
5000 for the planet. The disk is clearly seen to stand out in red, 
whereas the gap or the planet at 5 AU cannot be seen. The grid 
is written out in the beginning of the simulation and can be ex- 
plored in 3D using readily available open-source tools (e.g., par- 

Panel a) of Fig [ 



12 shows the integrated HCO + / = 7-6 emis- 



sion as predicted by LIME. This line has the frequency 624.2 
GHz and falls in the ALMA band 9. The model is, in this case, 
viewed face-on and placed in the distance of 100 pc correspond- 
ing to the typical distance to the closest disks. The resolution of 
the image is 0.005", which is about an order of magnitude higher 
than the best ALMA resolution. In the insert, we have zoomed in 
on the inner most 10 percent of the disk, where the protoplanet 
can be seen as a small dot inside the ring-shaped gap. The in- 
tensity variations across the innermost parts of the disk are due 
to the grid structure and not a property of the model. Such fluc- 
tuations are unavoidable when using random grids, but can be 
remedied by smoothing the image or adding more grid points in 
the affected regions. Panel b) shows the same image, but using 
the ALMA simulator (i.e., the simData task in the ALMA off- 
line data reduction package CASA) to generate simulated (u,v)- 
spacings and realistic noise. The present simulation is for a two 
hour exposure in a very extended ALMA configuration. Given 
that the mean distance between grid points on the surface is less 
than 0.01 AU, this simulation spans almost 5 orders of magni- 
tude in scales. 



In the panels c), d), and f) of Fig.[T2]we have flipped the disk 
so that we view it edge-on. In panel c) we show a model image 
for 183.3 GHz para-H20 / = 3i3-22o emission. This line falls in 
the ALMA band 5 and it is the only water transition available to 
ALMA. Only 6 band 5 receivers are planned so far, but here we 
show simulations of both a full array of band 5 receivers (panel 
d) and for the 6 planned receivers (panel f). In this edge-on view 
we see the lack of molecules in the mid-plane due to the simu- 
lated freeze-out. Both the LIME model and the ALMA simula- 
tions has been continuum subtracted so that the water emission 
stands out. In the case of a full array of band 5 receivers (panel 
d), the water emission is resolved spatially. This is not the case 
when using only 6 antennas for the simulation (panel f) because 
the antennas needs to be placed in a much more compact config- 
uration. 

This is a very highly opaque model, with the optical depths 
going up to several thousands. This is handled quite well by 
LIME although convergence takes 2-3 times more iterations than 
for HCO + . Because of the relative low abundance in this disk 
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Fig. 11. Visualization of the upper right quadrant of the disk model described in Sect. |4| Panel a) shows the H2 density with the 
insert showing the gap (shell) carved out around 5 AU. Panel b) shows the temperature and c) shows the molecular density which 
has a discontinuity around 90 K where water freezes out. The last panel, d), shows a cut through the grid, color coded according to 
the density. 



model, the water line does not maser in this particular example. 
Care should be taken when modeling potential maser lines with 
LIME because it does not handle masering accurately. LIME will 
produce a warning if the populations get inverted. 

The final example shows the multi-line ray-tracing capabil- 
ity of LIME. Here we calculate the spectrum of both flavors of 
methanol (type A and E) and SO2 around 338 GHz. We use the 
same disk model setup as above, but the lines are optically thin 
and we show the average spectrum over the disk so that result is 
relatively geometry independent. Figure 12 panel e) shows the 
resulting spectrum in a 2 GHz windows centered around 338 
GHz. In reality, methanol has many more lines in this window 
than what are shown here, but we are limited to the transitions 
for which the collision rates are know. 



5. Conclusions 

In this paper we have presented a new algorithm for solving 
the molecular excitation and radiation transfer problem in an ar- 
bitrary three dimensional geometry. The code uses a weighted 
stochastic point process to generate a random grid point distribu- 
tion and its corresponding Delaunay triangulation on which the 



photon transport takes place. The code is well suited for calcu- 
lating models with a large density contrast, which is particularly 
important for models of ALMA observations with a very high 
spatial resolution. Furthermore our code handles overlapping 
lines which is also going to be needed for proper modeling of 
ALMA data. Although emphasis has been put on making LIME 
capable of modeling ALMA data, the code is also well suited 
for modeling data from other (sub-) millimeter interferometers 
(e.g., SMA and CARMA) as well as data from single-dish ob- 
servatories (e.g., JCMT, APEX, Herschel Space Observatory). 

In this paper we show a comparison between the LIME code 
and a number of other molecular excitation and radiation transfer 
codes and we have applied LIME to a model of a protoplanetary 
disk with a protoplanetary condensation. 
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